Importation des données

En premier lieu, on importe les données, puis on assemble les deux fichiers en un seul jeu de données, dans lequel on modifie les noms de variables (étape \(1\) et \(2\) du TP).

pen.tra = read.table("pendigits.tra", sep = ",")
pen.tes = read.table("pendigits.tes", sep = ",")
pen = rbind(pen.tra, pen.tes)
names(pen) = c(paste(c("X", "Y"), rep(1:8, each = 2), sep = ""), "chiffre")
pen$chiffre = factor(pen$chiffre)

Un peu de dénombrement

Chaque chiffre est écrit en moyenne 1099 fois, avec quelques variations entre les chiffres

occur.chiffre = setNames(data.frame(table(pen$chiffre)), c("chiffre", "occurences"))
knitr::kable(occur.chiffre)
chiffre occurences
0 1143
1 1143
2 1144
3 1055
4 1144
5 1055
6 1056
7 1142
8 1055
9 1055 
ggplot(occur.chiffre, aes(chiffre, fill = chiffre)) + geom_bar(aes(weight=occurences)) +
    ylab("")

Premier exemple de chaque chiffre

Pour répondre à la question \(3\), il nous faut une fonction permettant de dessiner un chiffre.

dessineChiffre <- function(v, titre = NULL) {
    if (is.data.frame(v))
        v = unlist(v)
    don = data.frame(
        x = v[seq(1,15,by=2)],
        y = v[seq(2,16,by=2)],
        position = 1:8
    )
    g = ggplot(don, aes(x, y)) + xlim(0, 100) + ylim(0, 100) +
        geom_path(col = "gray50") +
        geom_text(aes(label = position)) + 
        theme_void() 
    if (!is.null(titre))
        g = g + ggtitle(titre)
    g
}
# pour tester : dessineChiffre(pen.tra[1,1:16], pen.tra[1,17])

Voici donc le premier exemple de chaque chiffre. On utilise ici le package gridExtra (et la fonction marrangeGrob()) pour pouvoir mettre plusieurs graphiques sur la même page. Et on se base sur le fait que la fonction ggplot() renvoie un objet stockable, et donc la fonction dessineChiffre() aussi.

premier = lapply(0:9, function(ch) {
    prem = subset(pen, chiffre == ch)[1,]
    dessineChiffre(prem, ch)
})
marrangeGrob(premier, nrow = 2, ncol = 5, top = "")

Tracé moyen de chaque chiffre

Dans un premier temps, nous devons calculer les coordonnées moyennes de chaque point, pour chaque chiffre (réponse à la question \(5\)).

pen.means = apply(pen[-17], 2, tapply, pen$chiffre, mean)
knitr::kable(round(pen.means))
X1 Y1 X2 Y2 X3 Y3 X4 Y4 X5 Y5 X6 Y6 X7 Y7 X8 Y8
0 35 86 12 58 15 20 51 7 86 31 89 68 59 89 22 75
1 15 61 44 78 70 90 77 80 68 54 48 33 45 16 60 1
2 18 77 42 99 67 80 51 46 20 19 12 9 53 5 99 4
3 25 84 57 100 87 85 65 61 82 43 91 17 50 2 3 6
4 43 100 22 79 6 51 43 40 85 50 86 60 71 31 63 0
5 41 91 43 76 57 59 36 29 26 33 38 50 43 58 59 60
6 88 99 52 87 21 58 7 27 33 3 81 11 62 31 11 23
7 3 91 45 98 79 81 71 47 53 15 34 18 40 34 81 34
8 57 82 40 80 52 52 51 24 35 17 40 37 68 68 49 81
9 69 81 53 83 45 81 57 83 79 71 90 43 61 14 18 5

Puis nous pouvons donc avoir la représentation des tracés moyens (question \(6\)), dans laquel nous voyons que presque tous les chiffres sont lisibles, et d’autres beaucoup moins (5 et 9 par exemple).

trace.moyen = lapply(0:9, function(ch) {
    dessineChiffre(pen.means[ch+1,], ch)
})
marrangeGrob(trace.moyen, nrow = 2, ncol = 5, top = "")

Comparaison des coordonnées pour chaque chiffre

Pour cela, nous pouvons représenter la distribution des coordoonées, pour chaque chiffre, avec des boîtes à moustaches par exemple.

distrib.coord = lapply(names(pen)[-17], function(var) {
    ggplot(pen, aes_string("chiffre", var, fill = "chiffre")) + geom_boxplot() +
        ggtitle(var) + 
        theme(legend.position = "none") + ylab("") + xlab("")
})
marrangeGrob(distrib.coord, nrow = 2, ncol = 2, top = "")

Représentation des tracés dans un espace restreint

Pour pouvoir avoir une vision synthétique (et répondre à la question \(8\)), il est possible de projeter les tracés dans le premier plan factoriel produit par l’ACP standard. On remarque que les deux premiers axes permettent d’obtenir 50 % de l’information.

acp = PCA(pen, graph = FALSE, quali.sup = 17)
ggplot(acp$eig, aes(factor(1:16), weight=eigenvalue)) + 
    geom_hline(yintercept = 1, linetype = 2, col = "gray50") + 
    geom_bar() + xlab("Composantes") + ylab("")

ggplot(setNames(acp$eig, c("eigen", "pct", "pctcum")), aes(1:16, pctcum)) + geom_line() +
    ylim(0, 100) + ylab("") + xlab("Composantes")

On peut donc représenter les tracés sur le plan factoriel.

ggplot(data.frame(acp$ind$coord, chiffre = pen$chiffre), aes(Dim.1, Dim.2, col = chiffre)) +
    geom_point()

Ici, il est difficile de comprendre les différences entre les points pour deux raisons :

  • il y a 10 chiffres et donc 10 couleurs, qui ne peuvent être assez différentes pour l’oeil humain
  • certains chiffres se ressemblent et donc les points se chevauchent

Une façon de résoudre ce problème (qui est le but de la question \(9\)), et de représenter les tracés de chaque chiffre séparemment.

ggplot(data.frame(acp$ind$coord, chiffre = pen$chiffre), aes(Dim.1, Dim.2, col =  chiffre)) +
    geom_point() + facet_wrap(~chiffre, nrow = 2) + 
    theme(legend.position = "none")

On voit ainsi que certains chiffres sont très éclatés (tel le 5, le 8 et le 0 par exemple).

Combien de tracés pour chaque chiffre ?

Finalement, dans la question \(10\), nous nous demandons si certains chiffres ont un tracé unique ou plusieurs tracés possibles. Pour cela, nous devons faire de la classification puis trouver un nombre de classes adapté. Pour cela, nous allons nous aider de trois choses :

  • Dendrogramme produit par la CAH, avec distance euclienne et critère de Ward
  • Evolution du \(r^2\) et du \(PseudoF\) sur les résultats de \(k\)-means, pour \(k=1,\ldots,10\)

Calculs des partitions

On va faire ces calculs pour chaque chiffre.

lk = 1:10
classif = lapply(0:9, function(ch) {
    sub = subset(pen, chiffre == ch)[-17]
    cah = hclust(dist(sub), "ward.D2")
    km = lapply(lk, function(k) {
        return(kmeans(sub, k))
    })
    I = km[[1]]$totss
    W = unlist(lapply(km, function(r) return(r$tot.withinss))) 
    B = unlist(lapply(km, function(r) return(r$betweenss))) 
    r2 = B / I
    PsF = (r2 / (lk - 1)) / ((1 - r2) / (nrow(sub) - lk))
    PsF[1] = NA
    return(list(cah = cah, km = km, r2 = r2, PsF = PsF))
})

Recherche du nombre de classes

Ce qui nous donne pour chaque chiffre les résultats suivants :

classif.plot = lapply(classif, function(res) {
    return(list(
        ggdendrogram(res$cah),
        qplot(lk, res$r2, geom = "line") + ggtitle("r2") +
            ylab("") + scale_x_continuous(breaks = lk),
        qplot(lk, res$PsF, geom = "line") + ggtitle("Pseudo-F") +
            ylab("") + scale_x_continuous(breaks = lk)
    ))
})
for (ch in 0:9) {
    print(marrangeGrob(classif.plot[[ch+1]], nrow = 3, ncol = 1, top = paste("Chiffre :", ch)))
}
## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

Proposition du nombre de classes

On peut proposer les nombres de classes suivant pour chaque chiffre :

k = c(2, 3, 3, 2, 2, 6, 2, 4, 5, 6)
knitr::kable(data.frame(chiffre = 0:9, "k choisi" = k, check.names = FALSE),
             align = c("c", "c"))
chiffre k choisi
0 2
1 3
2 3
3 2
4 2
5 6
6 2
7 4
8 5
9 6

Description des classes

Ainsi, voici les différentes possibilités de tracés pour chaque chiffre

for (ch in 0:9) {
    kk = k[ch+1]
    part = classif[[ch+1]]$km[[kk]]$cluster
    tmoy = data.frame(classif[[ch+1]]$km[[kk]]$centers, id = 1:kk)
    temp = data.frame(
            setNames(melt(tmoy[,c(seq(1,15,by=2),17)], id.vars = "id"), c("id", "xvar", "x")),
            setNames(melt(tmoy[,c(seq(2,16,by=2),17)], id.vars = "id"), c("id", "yvar", "y"))
        )
    sub = data.frame(acp$ind$coord[pen$chiffre == ch,], part = part)
    grid.arrange(
        ggplot(temp, aes(x, y, group = 1)) + geom_path() + facet_grid(.~id) ,
        ggplot(sub, aes(Dim.1, Dim.2)) + geom_point() + facet_grid(.~part)
    )
}